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ABSTRACT 


A two-dimensional model of a single channel unreflected circulating fuel reactor is 
used to investigate the effect of a radial temperature distribution on the neutron kinetics 
A slug flow model with no delayed neutrons is used. The effects of a nonuniform radial 
heat source distribution and radial heat transfer are taken into account. 



HIGH-POWER KINETICS OF A TWO-DIMENSIONAL CIRCULATING- 


LIQUID-METAL-FUEL REACTOR* 
by Michael J. Kolar 
Lewis Research Center 

SUMMARY 

A two-dimensional model of a single channel unreflected circulating fuel reactor is 
used to investigate the effect of a radial temperature distribution on the neutron kinetics. 

A slug flow model with constant fluid properties describes the heat transfer and fluid flow. 
A one-neutron velocity model with no delayed neutrons is used. The effects of a nonuni- 
form radial heat source distribution and radial heat transfer are taken into account. The 
results show that, in many practical cases, radial heat transfer has little effect on the 
kinetics. Furthermore, if changes in multiplication are not too large, the flux and radial 
temperature distribution maintain their steady-state shapes throughout the transient. For 
certain combinations of the physical parameters, resonances can occur in the time 
dependent flux and temperature. These resonances can be predicted by either a line 
reactor model or a two-dimensional reactor model; however, a two-dimensional model 
is necessary to predict the value of the steady- state flux and temperature following a 
transient. 


INTRODUCTION 

Since the early days of reactor technology, liquid -metal -fuel reactors have been the 
source of much speculation and research. Although this type of reactor was first sug- 
gested in 1941 (unpublished data by H. Halban and L. Kowarski mentioned in ref. 1), it 
received little attention until about 1947. At that time a program was started at Brook- 
haven National Laboratory to develop a Liquid-Metal-Fuel Reactor (LMFR). Basic 
research on the LMFR was continued until about 1957. During this time, a number of 
design studies were conducted (refs. 2 to 4) which indicated that this reactor concept was 
attractive economically. A number of LMFR experiments were planned in order to prove 
that the concept was technically feasible. The investigation was discontinued in 1957. 

^Information presented herein was offered as a thesis in partial fulfillment of the 
requirement for the degree of Doctor of Philosophy to Case Western Reserve University, 
Cleveland, Ohio, June, 1968. 


Since the LMFR is economically attractive, it is likely that it will be investigated 
again for use as a ground based power reactor. Furthermore, such a reactor could prove 
useful for space power generation since it could be operated as a relatively compact fast 
reactor. While there are no current programs to construct such reactors, there are a 
number of interesting research problems which can be done independent of such programs. 

Before discussing the particular problem of interest, it will be useful to describe 
two possible LMFR concepts. The first concept is that of a single fluid thermal reactor. 
The fissionable material is dissolved in a single fluid carrier (uranium in bismuth is one 
possibility) which circulates inside a closed loop. In the reactor portion of the loop, 
there is a block of graphite (or other) moderator material with cylindrical fuel passages 
in it. A neutron reflector surrounds this portion of the loop. The fuel solution flows 
through the reactor where it is heated and then through an external heat exchanger. The 
second LMFR concept is similar to the first except that no moderator or reflector is 
present. Such a reactor would operate in the intermediate or fast regions of the neutron 
spectrum. The reactor portion of the loop would be slightly larger in radius than the rest 
of the loop to reduce the neutron losses due to leakage. Nearly all other LMFR designs 
are variations of these two concepts. 

All previous research on the LMFR has dealt with the thermal reactor concept. In 
these studies, many authors have used an idealized model to study the reactor physics. 
This model consists of a single channel loop in which the cross-sectional area is con- 
stant throughout and the unreflected reactor is specified to be in some fixed portion of the 
loop. This model permits the study of certain aspects of the reactor physics without 
tying the problem to a specific system. Since it is the intent of this report to discuss a 
reactor physics problem which is common to both fast and thermal reactors, the single 
channel model will be employed. 

A reactor using a circulating fluid fuel has a number of advantages over one with 
stationary fuel elements. First, the structure of the reactor is relatively simple, that 
is, the fuel can be cooled in an external heat exchanger separate from the reactor core. 
Hence, the nuclear requirements in the core and the heat transfer requirements in the 
heat exchanger need not both be satisfied at the same location. This permits high spec- 
ific power design since, for example, materials such as tungsten (which has a high 
neutron capture cross section) can be used in the heat exchanger without affecting the 
nuclear requirements of the reactor. Fuel handling, fuel reprocessing, and waste dis- 
posal are easier in eirculating fuel reactors than in stationary fuel reactors. Another 
advantage of circulating fuel reactors is that fission products can be removed continu- 
ously. The removal of poisons (neutron absorbers) improves the neutron economy and 
permits higher fuel burnup. Consequently, less radioactive material is present in the 
reactor so that the potential hazard is decreased. One of the most advantageous features 
of a fluid fuel reactor is its inherent safety and ease of control. A liquid fuel which 
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expands on heating results in a negative temperature coefficient of reactivity. Since the 
rate of expansion is limited only by the speed of sound in the liquid, this effect is essen- 
tially instantaneous and tends to make the reactor self -regulating. Furthermore, a 
circulating fluid fuel reactor can be controlled by adjusting the fuel concentration at the 
reactor inlet. 

Circulating fluid fuel reactors also have several disadvantages relative to stationary 
solid fuel reactors. Changes in the fuel density or concentration can bring about fluctu- 
ations in the reactivity; this may occur, for example, if the local temperature changes 
in the fuel. Another disadvantage is that some delayed neutrons are lost from the core. 
This shortens the time between neutron generations and reduces the effectiveness of 
external control systems. Furthermore, there is considerable induced activity in pumps 
and heat exchangers. The induced activity produces radiation levels which require 
development of remote maintenance techniques. One of the biggest problems involved in 
all types of fluid fuel reactors is that of corrosion and erosion of structural materials by 
the fuel solution. 

Some of these disadvantages will be eliminated as experience is gained in handling 
liquid metals and in designing remote maintenance equipment. On the other hand, the 
problems of reactivity fluctuations and loss of delayed neutrons are inherent in circu- 
lating fuel reactors; thus it must be assumed that these reactors can be designed to 
operate efficiently and safely in spite of these disadvantages. A large portion of the liter- 
ature on circulating fuel reactors and the LMFR deals with these two problems. One 
may classify these as reactor kinetics problems since they involve the time -dependent 
behavior of the flux after some change occurs in the system. The reactor kinetics will 
now be considered in more detail. 

When the number of neutrons produced in a reactor changes from generation to 
generation, the neutron flux will change and a transient is said to occur. This situation 
is described by saying that a change in the multiplication or in the reactivity has occurred. 
After a change in multiplication in any reactor which has a negative temperature coef- 
ficient of reactivity, the neutron flux must attain a new steady state after some average 
temperature has been reached. It is often found, however, that the temperature rise and 
neutron flux can be out of phase so that flux and temperature oscillations can occur. 

Local temperature variations can have a large effect on the multiplication and, therefore, 
the spatial temperature distribution in the reactor must be known. 

In a circulating fuel reactor, the time history of the neutron flux after a change in 
multiplication depends not only on the temperature, but also on the fuel velocity distri- 
bution; this is true for two reasons. First, the velocity and temperature distributions 
are directly related through the energy balance in the reactor. Second, the fuel velocity 
affects the neutron kinetics because of the loss of delayed neutron emitters. If, for 
example, the reactor is operating at steady state with a certain fraction of the delayed 
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neutrons being lost to the system, then a change in the fuel velocity will change this frac- 
tion and a steady state will no longer exist. 

The fuel velocity in a circulating fuel reactor depends, in general, on the desired 
outlet temperature. For example, with a specified outlet temperature and a fixed inlet 
temperature, low velocities are necessary if the reactor power is low and high velocities 
if the power is high. While there is a complete spectrum of velocities which could be 
encountered in a circulating fuel reactor, there are many practical situations in which one 
of these two extremes (high or low power) exist. In each of these extremes, there are a 
number of simplifying circumstances which permit a rather clear exposition of the 
reactor kinetics. 

In a low power case, temperature changes in the system may be small. In this case, 
temperature effects on reactivity can be ignored and the kinetic effect of delayed neutrons 
can be studied alone. Fleck (refs. 5 and 6) has described a single channel circulating fuel 
reactor operating at low power by means of several linear partial differential equations. 
Expressions for the time dependent neutron flux and delayed neutron precursor distribu- 
tion are developed and the in-hour equation for the circulating fuel reactor is obtained. 

The fuel is assumed to move in slug flow, and only axial variations in the flux and pre- 
cursor distribution are considered. Expressions for the in-hour equation for a circu- 
lating fuel reactor have also been developed by Ergen (refs. 7 and 8) and Wolfe (ref. 9). 
The low power kinetics of circulating fuel reactors have been further discussed by 
MacPhee (ref. 10) and LoSurdo (ref. 11). The results of these studies show that the 
change in multiplication due to a change in fuel velocity is proportional to the square of 
the multiplication necessary to maintain the system at steady state, the delayed neutron 
fraction, and the fuel transit time across the reactor; it is inversely proportional to the 
velocity and the fuel transit time around the entire loop. 

When a circulating fuel reactor operates at high power, the fuel velocities encoun- 
tered are usually high enough so that a large fraction of the delayed neutrons is lost from 
the system. In this case it is often assumed that the reactor must operate on prompt 
neutrons alone. A reactor which is prompt critical has a period which is in the order of 
microseconds to milliseconds. A period this short may preclude the use of external con- 
trol systems; thus, the temperature coefficient of reactivity is the chief stabilizing factor 
in the system. A number of authors have studied the high-power kinetics of a circulating 
fuel reactor; in most cases a single channel model was used and all delayed neutrons 
were assumed lost to the reactor. Some of these analyses predict the time behavior of 
the flux and temperature as a function of axial position. A slug fluid flow model is 
assumed, along with constant fluid properties and constant flow rate. Ergen (refs. 12 
and 13), Weinberg (ref. 14), Welton (refs. 15 and 16), Thompson (ref. 17), and Nohel 
(ref. 18) have formulated the problem in terms of equations involving integrals over the 
past history of the neutron flux. This leads to a set of equations in terms of the time 
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variable alone. Fleck (ref. 19) has started with the partial differential equations for the 
neutron and energy balance. He has shown that, by expanding the flux in a time -dependent 
Fourier series, the resulting equations for the time -dependent coefficients are simply 
the relations used by the authors mentioned previously. The results of these analyses 
show that the negative temperature coefficient of reactivity is usually sufficient to stabilize 
the reactor after a step change in multiplication. However, if the reactor power, tem- 
perature coefficient, and the fuel transit time across the reactor are related in a certain 
manner, large oscillations of the flux and temperature can occur. 

The previous discussion indicates that a number of simplifying assumptions have 
been used in the study of the kinetics of circulating fuel reactors. Unfortunately, the 
conditions for these assumptions to be valid are neither quantitiatively defined nor uni- 
versally met. In past studies of the high power kinetics of the LMFR, the radial flux and 
temperature distributions have been assigned through a heuristic argument. Consider, 
for example, reference 19 where the radial flux is assumed to be exactly the same as the 
steady-state flux in a bare cylindrical reactor but the radial temperature distribution is 
assumed to be a constant. This will be referred to as a line reactor model. These 
assumptions are motivated by the physical situation which is often encountered in a cir- 
culating fuel reactor. From a mathematical point of view, these assumptions are incon- 
sistent since a radial flux distribution implies a radial heat source distribution; hence, 
the radial temperature distribution should not be a constant. Therefore, in the interest 
of obtaining a mathematically consistent solution and as a direct extension of the work of 
reference 19, a study of the high-power kinetics of a cylindrical circulating -fuel reactor 
has been made in which the radial terms in both the flux and temperature equations are 
retained; that is, a two-dimensional model has been used. 


SYMBOL LIST 


a length of reactor, cm 

a p-ja/h, dimensionless 

a z unit vector in z -direction 

2 2-2 
B geometric buckling, B^, cm 

B^ physical buckling, (t-S f p th g th - 2 a )/D, cm 

2 2 
B ^k geometric buckling of order n,k, cm 

C i constant 
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D 

E 

F 

F ii 

f 

f l> f 2 

G 

». 

G 

pi 

nzm 

g th 

h 2 

I 

Jn 


heat capacity of fuel solution, J/ (g)(°C) 
neutron diffusion coefficient, cm 
energy release per fission, J 
arbitrary function of r and z 

P l (B 2 - BjjJ/B 2 + 6kQ, dimensionless 

arbitrary function of time 
constants 

1 -a £ ?x o 

r M 111 (pc p )“ K 111 (v Q x)e u(r a - x), cm /J 


derivative of G with respect to t, cmy (J)(sec) 

/ T a -a x 

V 1 


K nlm^ e 1 sec 


fast nonleakage probability 
P l /B 2 , cm 2 

/ OO 

S(x) sin (J2x) dx, dimensionless 

V 

zero-order Bessel function 
first-order Bessel function 


^nZm 

I/c K 

‘ < z > z m < z - v, 

K 1 

vQo“m. 

cm 2 

K 2 

T le X /< 3 - 

Ti ex ), dimensionless 

k 

finite medium multiplication 

k 

OO 

infinite medium multiplication 


dimensionless 


L diffusion length, cm 

l ^ finite medium lifetime for neutrons with speed v k , sec 
l q infinite medium lifetime for neutrons with speed v k , sec 




N 

integer > 3 

P 

reactor power density, J/ (cm )(sec) 

• 

P 

3 2 

derivative of P with respect to t, J/ (cm )(sec ) 

P L 

thermal leakage probability, 1 - P^ 

P NL 

2 2 

thermal nonleakage probability, (1 + L B ) 

P 0 

initial reactor power density, ESj$q, J/(cm )(sec) 

P 

Laplace variable, sec - * 

Pth 

resonance escape probability 

2 

Pi 

6* 0 - A* 

Q 

function of time, In (P/Pq), dimensionless 

Q 0 

ES f (pc ) _1 , (cm 2 )(°C) 

Q 

-2 

second derivative of Q with respect to t, sec" 

R 

reactor radius, cm 

^(r) 

J o(£j r ), dimensionless 

r 

radial coordinate, cm 

S 

-G(t)/G(0), sec -1 

T 0 

inlet temperature, °C 

T 1 

2 

y l u / Pi> dimensionless 

T lex 

T, when u = u , dimensionless 
x ex 

t 

time, sec 

t l» t 2 

arbitrary time, sec 


ft) where t < x 

u(t - x) 

unit step function, ^ where t>x 

v k 

neutron speed, cm/ sec 

v O 

velocity of fuel solution, cm/ sec 

X 

dummy integration variable 

y 

du/dz, °C/cm 
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Z n< z > 

z 

z 

a 


0 

7 

?! 

6(t-x) 


5 ii 


6k 

6k 


0 


h 

e 

0 line 


e 

e 

e 

e 


max 

R 

0 


2-d 


K 


V 


V i 

k 

P 

S„ 


■a 


sin (e n z) 

axial coordinate, cm 

Pjz/h, dimensionless 

thermal diffusivity, n(pc)~ , cm /sec 

r 


2 1 

4>l/ dimensionless 


negative of temperature coefficient of multiplication, C 
r M m , °C' 1 


O/-,- 1 


Dirac delta function defined by 



f(x) 6(t-x) dx where t^<t < t 2 


JO where i # j 

Kronecker delta, 1 , , „ 

’ (1 where i = j 

change in multiplication, k - 1 

initial change in multiplication 

njr/a (where n = 1, 2, . . . ), cm' 1 

Jl(^l)/Jji(^ 2 )> dimensionless 

temperature rise above inlet temperature, °C 

equivalent to 0 Q but developed from line reactor model of circulating fuel 
reactor, °C 

maximum temperature reached during excursion, °C 
wall temperature rise, °C 

initial and/or steady-state temperature rise, °C 
equivalent to 0 q, °C 

sum of molecular and turbulent conductivities of fuel solution, J/(cm)(sec)(°C) 

number of neutrons produced per fission 

solution of Jq(ia) = 0 

radial buckling, y./R, cm" 1 

density of fuel solution, g/cm^ 

macroscopic absorption cross section, cm" 1 

macroscopic fission cross section, cm" 1 

transit time for fuel across reactor, sec 
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ex 

*(t) 

m 

*( 0 ) 


<i> 


max 




0 

-k 

*1 

•k 

h 


1, l-term 
1 

1, 2 -term 


<P 
<P 0 

cp s (z) 

* 

<2 2 


function of axial coordinate, °C 
u evaluated at z = a, °C 

flux amplitude equivalent to <i>j, [(cm 2 )(sec)] - * 
derivative of <3>(t) with r espect to t, |(cm)(sec)] -2 
initial flux amplitude, [(cm 2 ) (sec)] - ^ 

maximum flux amplitude during excursion, |(cm 2 )(sec)]”* 
steady-state flux amplitude, [(cm^)(sec)] - ^ 


l , k coefficient of flux in Fourier-Bessel expansion, [(cm^)(sec)]'’*' 
derivative of w ith respect to t, [(cm)(sec)] _2 

solution for <f>j when one-term approximation is used, [(cm 2 )(sec)]"'*' 

solution for <3>* when one -axial -term, two-radial -term approximation is 
used, j(cm 2 )(sec)]‘^ 
neutron flux, |(cm^)(secj]'^ 
initial neutron flux, [(cm 2 ) (sec)] * 
steady state axial flux, [(cm 2 )(sec)] -1 
T a «^ 2 - ^j/2n, dimensionless 

G(0)Pq/?^j i , sec 


PHYSICAL MODEL 

The reactor consists of a right circular cylinder with radius R and length a with 
fuel moving from left to right as shown in figure 1. The radial coordinate is r; the axial 
coordinate is z, and symmetry about z is assumed. 

The fuel enters the reactor in a uniform thermodynamic state, that is, the external 
heat exchanger is assumed capable of maintaining a constant temperature at the reactor 
inlet. The temperature of the walls is assumed, in general, to be a specified function of 
the axial coordinate and time. If the flow is turbulent, the turbulent conductivity is taken 
as constant; in this case, time averaged quantities are used in the energy equation 
(ref. 20). Radial molecular conduction and turbulent heat transfer are therefore repre- 
sented by the same term in the energy equation. An internal heat source is present in 
the fuel because of fission heating, the heat source being taken as proportional to the 
neutron flux. 
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Figure 1. - Geometrical model of circulating fuel reactor; fuel 
temperature is Tg at z = 0. 


The reactor is an entrance region as far as the fluid flow is concerned. Because of 
the low Prandtl number of liquid metals, the fuel velocity can be considered constant 
over a cross section of the reactor. Thus, the assumption of slug flow seems reasonable. 
In the present discussion, the effects of free convection cells, viscous dissipation and 
axial conduction as compared to axial convection will be ignored. It is assumed that tur- 
bulence does not affect the neutron kinetics except through the energy equation. Also, all 
physical parameters involved in either the heat transfer or fluid flow are assumed to 
remain constant throughout the problem. 

A one-neutron velocity model is employed. It is assumed that all delayed neutron 
emitters are swept out of the reactor before decaying so that the reactor operates on 
prompt neutrons alone. The neutron flux is taken to be zero at the entrance, the exit, 
and the walls of the reactor. While changes in density are assumed negligible when 
considering the heat transfer and fluid flow problem, it is assumed that changes in den- 
sity can cause significant changes in the multiplication. These changes occur simultan- 
eously with changes in temperature. 

The reactor is initially operating at some power when a change in the multiplication 
is made; the change is uniform and instantaneous at all points in the reactor. The 
behavior of the neutron flux and fuel temperature after this change has been made will be 
studied. 


ANALYSIS 

Basis Equations 

The time dependent diffusion equation for one -velocity neutrons can be written 
(ref. 21) 
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( 1 ) 


-DV 2 <p(r,t) = ^S f P th g th <p(r,t) - i- - S a <p(r,t) 

k 

where cp( r, t) is the neutron flux and all other symbols are defined in the list of symbols. 
Since symmetry about z has been assumed, the flux is a function of only r, z, and t. 
The flux is zero at the reactor extremities, symmetric about the axis, and is a specified 
function of position at t = 0. Thus 


<p(r, 0,t) = 0 

(2a) 

<?(r, 

, a, t) = 0 

(2b) 

<p(R,z,t) = 0 

(2c) 

d<£_ 

= 0 

(2d) 

dr 

r=0 


(p( r,z, 0) = (p 0 ( r,z) 

(2e) 


Dividing equation (1) by and defining 
-1 a 
vS f p th S a result in 




and k = 

oO 


1 ° -j: _ fc^th " + L ' 


i 

r 



(3) 


H P NL s + l2b2 ) P L 5 1 - P NL’ k = k ooS th P NL’ and l th s 1 0 P NL’ ec * uation ( 3 ) 
becomes 


ifh - (k - + P 


at 


NL' 


Sr\ Sr / dz 2 


(4) 


where B 2 = (v-^/R) 2 +' (7r/a) 2 . Let k, the multiplication, be written as k = 1 + 5k where 
6k represents the deviation from equilibrium. Then the time-dependent diffusion equation 
finally becomes 
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l ^ = 6k <p + P T B -2 B 2 <p + I _1 (r l2\ + (5) 

th 3t L r ar V 9*7 9z 2_ 

All of the parameters which multiply the flux and its derivatives in equation (5) vary 
with temperature. However, it is assumed that only variations in 6k have a significant 
effect on the reactor kinetics. To understand the conditions under which this assumption 
is valid, it is necessary to consider one solution of equation (5). Assume all parameters 
in the equation are constant; then the solution is cp( r, z, t) = F(r, z)e^7th where F(r, z) 
is some function of the spatial coordinates. F(r, z) arises directly from the terms in 
brackets on the right side of equation (5). Note that only the ratio 6k/ 1 ^ is involved in 
the transient part of the solution for the flux. Next the parameters in equation (5) are 
permitted to vary slightly with temperature. The function F(r,z) is assumed to change 
very little under these circumstances relative to the change in the exponential. Hence, it 
seems reasonable to assume that temperature changes in those parameters involved in 
the transient solution for the flux will be more important than temperature changes in pa- 
rameters involved in the spatial part of the solution. If this is true, only small quantities 
of second order are neglected by assuming that the temperature enters equation (5) only 
through 6k. While the functional dependence of 6k on temperature is ordinarily quite 
complicated, it is usually possible to use a linear variation of 6k with temperature as an 
adequate representation. Hence, 6k is assumed to have the form (see ref. 21, pp. 308 
to 327) 

6k = 5k Q - y0(r,z,t) (6) 

where 6kp is an initial change in multiplication made instantaneously and homogeneously 
throughout the reactor and y, the negative of the temperature coefficient, is assumed 
constant throughout the transient (y is defined by the relation y - - 9k/ 9 9). Note that 
0(r,z,t) is the temperature rise above inlet temperature at the spatial location and time 
at which the change in multiplication occurs. 

Using equation (6) in equation (5) results in one equation which relates the flux and the 
local temperature rise 

_ dz _ 

A second expression relating <p and 6 can be obtained from the Navier-Stokes equa- 
tions. For an incompressible fluid in slug flow, the continuity and momentum equations 
become p = constant and v = v^a , respectively. The energy equation becomes (ref. 20) 
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( 8 ) 


pc (— + v 0 — ' Wr,z,t) = EE f <p(r,z,t) + k- — (r — \ 
p \St u dzj 1 r 9r \ 9r/ 

where axial conduction has been ignored. Equation (8) is simply an energy balance; 

E Y,^(p represents the fission heating, and k is the sum of the molecular and turbulent 
(if any) conductivities. In the case of turbulent flow, the temperature and heat transfer 
parameters are assumed to be time averaged quantities. 

The temperature is measured with respect to the inlet temperature which is constant. 
The temperature rise is symmetric about the z-axis, a specified function of axial posi- 
tion and time at the wall and is a specified function of position at t = 0. Therefore, 


0(r, 0,t) = 0 

(9a) 

30 

= 0 

(9b) 

3r 

r=0 


0(R, z,t) = e R (z,t) 

(9c) 

0(r,z, 0) = 0 Q (r,z) 

(9d) 


Define Qq = ES f (pc p ) 1 and a = /c(pCp)" 1 ; then equation (8) may be written 


(— + v o —) e = + 01 ~ — ( r — ) 

Vat U 3z/ u r ar V 0r/ 


( 10 ) 


Equations (7) and (10) are a coupled set of partial differential equations; together 
they describe the kinetics of a LMFR which behaves according to the assumed physical 
model. 


Fourier-Bessel Form of the Basic Equations 

Equation (7) is a second-order nonlinear partial differential equation, the nonlinearity 
being a result of the temperature coupling. The solution to such an equation depends on 
the magnitude of the nonlinear term when compared to the other terms in the equation. 

If equation (7) is used to describe a high power LMFR in which a reasonably large initial 
change in reactivity has occurred, then the equation is highly nonlinear. On the other 
hand, for a low power LMFR or for small perturbations about equilibrium in the high 
power case, the temperature rise may be small and the nonlinearity is correspondingly 
small. 

By a qualitative discussion of equation (7), a great deal of insight into the kinetic 
behavior of a reactor can be obtained. For example, if the temperature rise is taken to 
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be zero, the equation is linear and the neutron flux is given by a product of two space 
functions and an exponential in time. The period of such a system is proportional to 
l/6k (ref. 21). If the temperature is permitted to rise in time, however, it can be seen 
that eventually the initial change in multiplication will be exactly cancelled by the temper- 
ature term. At this point, the reactor period is infinite or, in other words, the flux no 
longer rises. This would occur if the temperature were coupled to the flux only through 
equation (7). 

Actually, the flux and temperature are also coupled through equation (10). If equa- 
tion (10) were solved for the flux and the result inserted in equation (7), a single equation 
for the temperature would be obtained. One of the terms in this equation would involve 
the second derivative of the temperature with respect to time. Hence, one might suspect 
that oscillatory solutions for the temperature are possible. If this were true, then flux 
oscillations would also occur. Any study of the reactor kinetics must therefore deter- 
mine whether such oscillations are possible and, if they are possible, whether the oscil- 
lations are bounded and are damped. It has been previously shown (ref. 19) that, for a 
line reactor model, such oscillations can exist and are bounded and damped. It is the 
purpose of the present work to investigate the effect of radial temperature gradients and 
heat transfer on the behavior of these oscillations. To do this, the coupled equations (7) 
and (10) must be solved. 

Two possible methods present themselves for solving these equations. A direct 
numerical solution by finite difference techniques is one possibility; this method would 
give little insight into the physics of the problem since each reactor system would be a 
special case. The second method and the one which will be used here is to solve equa- 
tion (10) for 0(r, z, t) by analytical methods and use the resulting expression in equa- 
tion (7). The flux is then expanded in a time -dependent Fourier -Bessel series and the 
important parameters are determined. 

An expression for 6( r, z,t) in terms of cp( r, z, t) can be obtained in the following 
manner. It is shown in appendix A that the application of a finite Hankel transform 
(ref. 22) with respect to r, followed by a Laplace transform (ref. 22) with respect to t 
to equation (10) results in an ordinary first-order differential equation for the transformed 
temperature. After solving the differential equation, the temperature is retransformed 
with the result 


0(r, z ,t) = — 

R 2 




J 0 (^j) 


( 11 ) 


where 
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e(%, z >t) = Q 0 


/ z/v 0 


(pit 


-a^x 

ij, z-VqX, t-x)e u(t— x) dx 


/ z/v 0 


-a£x 


9 o(£i> z— v Qx) e 6(t-x) u(t-x) dx 




/ z/v Q 


-a|jX 

0r( z -V(jX, t-x)e u(t-x) dx 


( 12 ) 


Here, <p(£j,z, t) is the finite Hankel transform of the flux cp, 0 q(|j, z) is the finite Hankel 
transform of the initial temperature distribution, and 0 R (z,t) is the specified wall tem- 
perature. Equation (11) can now be inserted into equation (7) to yield a single equation 
involving the flux. 

Assume that the flux can be expanded in a time dependent Fourier -Bessel series of 
the form 


<p(r, M) = 2 £ 4(t)Z ( Z )(ft.(r) (13) 

n=l j=l J 

where <i>^(t) is a time -dependent flux amplitude; Z n (z) = sin ( e n z ), e n = nw/a; and 
(R. (r) = J„(|.r), £. = vJR, v . is defined by the equation Jq(za) = 0. The form of the 
series expansion for the flux is suggested by the steady -state solution for the flux in a 
cylindrical reactor with no temperature coupling, that is, cp( r, z, 0) = $(0)Z^(z)(Rj(r) 

(see appendix D). 

It is shown in appendix B that, when equations (11) and (13) are used in equation (7), 
then the usual methods of Fourier analysis lead to the following coupled set of equations 
for the coefficients <f^(t): 
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I Ill 


= K + p l b_2(b2 - <&>] *?<*> 


-V ) M yk * ) n (t) I dxe 


n,j 

m, i 


2a ^ 


i 


-a^fx 


Qo^(t-x) K nZm ( v oX) 


*V 


a 

- 0 R (z-VoX,t-x)Z^(z)Z n (z) dz 


r 



-a^t , 

2e x ^(t) 


M ijk 2 


2-^ J (' t Z;(z)Z n (z)e 0 (ii,Z-T 0 t)dZ 0==t<T a 


n,J 
m. i 


~ ¥ o' 


(14) 


7 

Z th^ 


(t) = [ak 0 


+ P l B" 2 (B 2 - B 2 ^ $*(t) 



n ,j 

m, i 


“iJk^nW 




Q 0 «L^m^) 


20 !^. 


V 


,t-x)Z l (z)Z n (z) dz 


t > T a (!5) 


Here, 


K nZm< v O x )= 7 - 

C l 


f Z n^ Z l ^ z ) Z m^ z_v O x ^ dz 

*V 


and 
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TABLE I. - FIRST FEW VALUES OF K nim (v,jX) AND M yk 


n 

l 

m 




1 

1 

1 

A/l + icosH + AcosM 





77 \ 3 

a 3 a / 


1 

2 

1 

-1 f-Acos H - v /l - *\ sin “ + A cos. 

3T7x“] 




277 L 4 

a \ a/ a 4 

a J 

1 

1 

2 

— [1 - cos ^2 + 7 r /l - sin AEZEI 





277 L 

a \ a/ a J 


2 

1 

1 

J. f- A COS « - * /l - *\ sin 2 + A cos 

3ttx1 




277 L 4 

a \ a/ a 4 

a J 

1 

2 

2 

-1/i cos 

H + I! cos 2 I x + 4 cos 3 ttx\ 





277 \3 

a 15 a 5 a / 


2 

1 

2 

A/icos 

!E + 32 cos 2ttx + 4 CQS 3ttx\ 





277 \3 

a 15 a 5 a / 


2 

2 

1 

A (i + A® 

cos — + -A- cos 





77 \ 15 

a 15 a / 


2 

2 

2 

A/i-l 

cos2EL + AcoslH\ 





277 \ 3 

a 3 a / 






r & 


where 

m 

( V 0 X ) = 7- J 

j Z n ( Z ) Z 7 ( z ) Z m( Z ' V fl X ) dz 












* 


i 

j 

k 


M ijk 


1 

1 

1 


0. 37550 


1 

2 

1 


. 06953 


1 

1 

2 


-. 10606 


2 

1 

1 


. 06953 


1 

2 

2 


-.21405 


2 

1 

2 


-.21405 


2 

2 

1 


. 14225 


2 

2 

2 


-.04393 






r 


where 

M ijk : 

2 

Of - — 1 q 

- I r(R i (r)<Rj(r)(Rj j: (r) dr 





R [ J l("k)f 

Jo 



* /~1 

Values of / xJ 0 (^x)J 0 (^x)J 0 (y k x) dx were obtained 

from ref. 23. 
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Ft 

M ijk ° 2r 2 / r «i(r)« jfr^r) Or 

Table I lists the values of the first few m and M.j^. Equations (15) and (16) are an 
infinite coupled set of nonlinear integro-differential equations. It will be shown that in 
many cases of practical interest this set can be approximated by a single term. 

VALIDITY OF ONE-TERM APPROXIMATION 

The assumption is now made that 0 R (z,t) is zero for all z and t. With this as- 
sumption in mind, an argument is given in appendix C to the effect that, under certain 
conditions, the first term in the set (eq. (14) and (15)) will represent the kinetics problem 
extremely well. The basis of the argument is that the steady -state flux and temperature 
distributions obtained by taking either a single term in the set, two terms, or the exact 
steady-state solution when the conductivity is zero all yield the same result. This 
implies that the second term in the series for the flux is an extremely small perturbation 
on the first term and all other terms are smaller than the second. * It should be noted 
that a one-term approximation implies that the spatial flux distribution and the radial 
temperature distribution do not change during the transient. 

2 -2 

The conditions for a one term approximation to be valid are Skg « P^CjB , 

Sk Q « |P L (£ 2 - £ 2 )B~ 2 |, and r a « 1. Then the kinetic equations reduce to 

f ~ a Zi x 

Z th $(t) = 6k 0 $(t) - y*( t)M m 1 dx e Q 0 K in (v () x)$(t-x) 

•A) 



"''Note that the results obtained in appendix C make this statement plausible. In order 
to further investigate this statement, a numerical solution to the one -term equations was 
compared to a numerical solution of the exact equations given in ref. 19 for a line reac- 
tor. The two methods were very nearly the same in their results. 
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(17) 


l to*® = ak 0 *(t) - yM in $(t)Q 0 



-a^x 

dx e 1 $(t-x)K 111 (v Q x) 


t>T 


where the subscripts on #(t) have been dropped and the fact that ? 0 (£., z) contains the 
term 6^ has been anticipated. This is demonstrated in appendix D where the initial 
steady- state temperature distribution is derived. 

The kinetic behavior of the temperature is described by the equations 


6>(r,z,t) = 


R' 


'[fiO'ijp 


® t)J 0 (r 


(18a) 


0(1^ z,t) = Q 0 


t _ 

I (ptepZ-vtft, t-x)e A dx 


0(lpM) = Q 0 


f /v °_ 


+ B 0 te lf z-v 0 t)e 


-a^x 

z-VgX, t-x)e dx 


-a I? * 


0 < t < T a , t < z/v n (18b) 


0 < t < r a , t > z/vq, or t > 


(18c) 


STABILITY OF SYSTEM 

A reactor system is said to be stable if a departure from an equilibrium state ulti- 
mately results in a return to equilibrium. Linear stability implies that the disturbance 
which causes the departure from equilibrium is arbitrarily small, while nonlinear 
stability implies that the departure from equilibrium is large enough for nonlinear effects 
to become important. 

Consider the kinetic equation for the flux after a time greater than r (eq. 17): this 
equation describes the flux after a long period of time so that its stability should guaran- 
tee the stability of the system. Following Welton (ref. 16), define P(t) = E2,$(t) and 

2 * 

-1 -a£jX 

G(x) s yM in (pc p ) Km(VgX)e u(r a - x); then equation (17) can be written as 
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(19) 


J th P = pj6k 0 - G(x)P(t-x) dx 

Using the definitions P(t) = and S(x) = - G(x)/G(0) in equation (19) and diff er - 

entiating the result with respect to time give 

Q + fi 2 (e Q - 1) = fi 2 f°° dx S(x)[e Q(t_x) - l| (20) 

where fi 2 s G(0 )Pq/Z^.j 1 and the fact that S(x) dx = 1 is used. When the right side 

of equation (20) is zero, the resulting equation is that of a one -dimensional, nonlinear 

oscillator with no driving force or damping. The motion of such a system is strictly 

periodic. As a matter of fact, for small Q (i. e. , for linear power changes) and with the 

• • 2 

right side of equation (20) equal to zero, the result is Q + fi Q = 0 which shows explicitly 
2 

that fi plays the role of a frequency squared. This equation also shows that, for linear 
deviations from a steady state, it may be possible to obtain purely sinusoidal motion of 
the system. 

The right side of equation (20) acts like a hysteresis effect and, under suitable con- 
ditions, will lead to damping. It should be noted that if delayed neutrons are present, an 
additional term appears on the right side of equation (20) which always leads to damping 
(ref. 13). Welton has shown that, even in the absence of delayed neutrons, equation (20) 
is never unstable if 

~ S(x) sin (fix) dx > 0 (21) 

o 

for all possible values of fi. It is shown in appendix E that, if a£jT a is small enough 

-ol£t 

to permit the approximation e =* 1, then S(x) fulfills the condition (eq. (21)) pro- 

vided that fi < 3 tt/t . For larger values of fi, the left side of equation (21) oscillates 
between positive and negative values. Hence, it is concluded that the system is damped 
at least for values of fi < 3ir/ r_ . 

<x 

2 The procedure followed by Welton is to obtain the Lagrangian corresponding to the 
motion described by equation (20) from which the Hamiltonian is obtained. Now the 
Hamiltonian corresponds to the total energy excess above equilibrium. If the Hamiltonian 
is positive definite and bounded, then so is the energy thus implying conservative motion, 
and unconditional stability is therefore guaranteed. For many reactor systems this con- 
dition is actually too strong. 
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Fleck (ref. 19) has argued that equation (19) is bounded and is certainly damped in 

the case of nonlinear disturbances but that purely oscillatory motion can occur in the case 

of linear disturbances for certain values of fi. 2 Specifically, those values of fi which 

make the integral f °° G(x) cos (fix) dx equal to zero result in the flux and temperature 
*'0 

being completely out of phase. It is shown in appendix E that this is the same integral 

used in Welton's stability criterion. For the particular G(x) used in this report, the 

"resonance" values of fi, that is, those values of fi which make /'°° G(x) cos (fix) dx 

•'0 

equal to zero, are given by fi=N7r/T a where N = 3, 4, 5, . . . In most practical cases, 
it will be found that fi < 3jr/ r a ; when the equality holds, the fuel temperature rise will 
be so large that many of the previous assumptions about various parameters being con- 
stant will no longer be valid. Using the definition of fi and the fact that the steady -state 
flux amplitude is given by 4 >q = ff6kQ(yQ Q M.|^jT a )~l (see appendix C, eq. (C-10)) shows 
that resonances occur when 


T a 6k 0 _ 3N 2 tt 2 


th 


N = 3, 4, 5, . . . 


( 22 ) 


It is concluded that the system will be stable for the case of high -power kinetics for 
all disturbances which are nonlinear. However, damping may be quite slow near one of 
the resonances described by equation (22). It is interesting to note that equation (22) is 
identical for either the line reactor or a two-dimensional reactor model and is independ- 
ent of the radial distribution of internal heat sources. 


NUMERICAL RESULTS 


In order to compare the present work with the line reactor model and to verify the 
analytical results, a few numerical examples were run on a digital computer using the 
computer program Mimic (ref. 24). The one-term equations (eqs. (16), (17), and (18)) 
were used along with the parameters given in table II; these parameters are similar to 
those given in reference 19 for a Brookhaven National Laboratory U-Bi LMFR design. 
With the density and heat capacity held constant, the thermal conductivity was varied 


°Fleck has assumed that the boundedness and nonperiodicity of solutions to equa- 
tion (19) implies that the solution damps. By expanding the flux in a Fourier series and 
examining the individual components he is able to show that periodic motion is possible 
only for linear disturbances of the system and that these motions occur when 
/"°° G(x) cos (fix) dx = 0. 

Jo 
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TABLE H. - PARAMETERS USED IN NUMERICAL CALCULATIONS 


Length of reactor, a, cm 

170. 031 

Reactor radius, R, cm 

92 

Negative of temperature coefficient of reactivity, y, °C“ 1 

5X10' 5 

Diffusion length squared, L 2 , cm 2 

229 

Geometric buckling, B 2 , cm" 2 

1. 026X10' 3 

Sum of molecular and turbulent conductivities of fuel solution, k . 

>0, <0.155 

J/ (cm)(sec)(°C) 


Heat capacity of fuel solution, c p , J/(g)(°C)) 

0. 149 

Density of fuel solution, p, g/cm 3 

10 

Thermal diffusivity, a, cm V sec 

>0, <0. 104 

Transit time for fuel across reactor, r a , sec 

1 . 0 

Finite medium lifetime for neutrons with speed v k , £ th , sec 

6. 56X10" 4 

Initial flux amplitude, $(0), [(cm 2 )(sec)]” 1 

1.0X10 15 



Figure 2. - Comparison of line reactor and two -d ime ns ion ai reactor model 
with a = 0. 


between zero and 0. 155 joule per centimeter per second per °C to study its effect on the 
kinetic equations. The one -term equations for a line reactor were also programmed 
and the results compared to the two-dimensional reactor. 

Figure 2 shows a plot of the centerline outlet temperature as a function of time for 
the line reactor and for the two-dimensional reactor with the thermal conductivity equal 
to zero. The figure indicates that the temperature history behaves nearly the same in 
each case but that the absolute values of the temperature rises are quite different. To 
see why the difference in magnitude exists, the analytic expressions for the steady -state 
outlet temperature and the steady-state flux amplitudes for the line reactor and the two- 
dimensional reactor can be compared. These are 
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i am 


inpaiaiiii m 





7 line (0 ’ a ' Z Z~’ 

ir ' 


5k n jr 

($ ° )line 


0 a 


and when a = 0 


0 2-d(°’ a ) " 


26 k Q J 1 (i/ 1 ) 

* M 111 


-d 


6k () ir J i(i/ 1 ) 

r<Va M m 


It is seen that the radial distribution of heat sources, which gives rise to the term, 

causes the difference in magnitude in temperature. The same result is obviously true for 
the flux. It is concluded, therefore, that a change in the radial internal heat source dis- 
tribution will alter the magnitude of die temperature in the reactor but it will not disturb 
the kinetic behavior of the reactor significantly. 

Figures 3, 4, and 5 are plots of the flux amplitude and axial outlet temperature for 
several values of Skg. These correspond to "values of N" of 0. 632, 1.28, and 3.22 
where N is computed from equation (22). These "values of N" indicate how close the 
system is to one of the predicted resonances. In figure 3 the system is practically in 
phase and the flux and temperature simply approach their steady -state values along a 
smooth curve. In figure 4 the system is closer to its first resonance and the flux and 
temperature both oscillate about their steady -state value before damping. Note that the 
temperature oscillations are less severe than the flux oscillations. This effect is more 
obvious in figure 5 where the system is very near its first resonance. 

The condition for a resonance was obtained by ignoring the effects of thermal con- 
ductivity. Figures 3, 4, and 5 were obtained from the two-dimensional reactor kinetics 



Figure 3. - Flux amplitude and centerline outlet temperature for "N" = 0.632; 
6k 0 = 0.001. 
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Kt-s-fi nl feslsl i i 1 ■ 




Figure 4. - Flux amplitude and centerline outlet temperature for "N" = 1. 28; 
6k 0 = 0.004. 


$/$max where ^max = 8 *93xl0 17 [(cm 2 )(sec)] 1 

8/e max where 6 max - 1.668xl0 3 ° C 



Figure 5. - Flux amplitude and centerline outlet temperature for "N" - 3. 22; 
6k 0 = 0. 025. 


equations with k = 0. 155 joule per centimeter per second per °C. It is concluded, 
therefore, that the position of the resonances of the system is not affected by radial heat 

9 

transfer provided, at least, that 1 t « 1. Three runs were made with the conditions 
used in figure 4 while the thermal conductivity was given values of 0. 0, 1. 55x10”°, and 
1.55x10” joule per centimeter per second per °C. (This is equivalent to holding k 
constant and varying R in the term a£*r.) The results showed that the kinetics were 
nearly identical, the only change being a difference of 0. 002 C out of about 222° C in the 
outlet temperature between the lowest and highest values of k. Hence, it is expected 
that, in many cases of practical interest, radial heat transfer can be ignored. 
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SUMMARY OF RESULTS 


The following conclusions can be drawn about the high-power kinetics of a single 
channel unreflected circulating fuel reactor with a nonuniform radial heat source distri- 
bution and radial heat transfer permitted. These conclusions are valid provided the one- 
neutron velocity diffusion equation is a valid description of the system and a slug fluid flow 
model is adequate. 

If a step change in multiplication is made homogeneously throughout the reactor and 
is such that 6kQ « P L ejB" 2 , 6kQ«P L |£ 2 - £||b~ 2 , and r a « 1 where 6kQ is 
the change in multiplication, P^ is the thermal leakage probability, B 2 is the buckling, 
and ^ are first and second modes of the radial buckling, a is the thermal dif- 
fusivity, and r a is the transit time across the reactor, then the following hold: 

1. The neutron flux will essentially maintain its initial steady -state radial distri- - 
tion throughout the ensuing transient. 

2. The temperature will essentially maintain its initial steady-state radial distri- 
bution throughout the transient, provided the wall temperature rise is zero. 

3. A number of oscillations in the flux and temperature are possible; however, a new 
steady state will be reached. These oscillations will be bounded but will have a higher 
frequency and will be weakly damped when a certain relation exists between the param- 
eters of the problem, that is, TgSkgZ^ = 3N v /8 where N = 3, 4, 5, . . . and is 
the neutron lifetime. 

4. A line reactor model gives an adequate description of the flux and temperature 
time histories, but the two-dimensional model must be used to obtain the magnitude of 
the flux and temperature. 

5. In many cases of practical interest, radial heat transfer can be ignored. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, July 3, 1968, 

122-29-03-01-22. 
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APPENDIX A 


SOLUTION OF ENERGY EQUATION IN TERMS OF THE FLUX 


Consider the energy equation in the following form: 


(— + ^ 0 = Qq<p + a — — (r — 

vat u 3z / u r 3r \ Sr, 


de'j 


( 10 ) 


Apply a finite Hankel transform (ref. 22) to this equation letting 


/*R 

0(%,z,t)s re(r,z,t)J 0 (r^.) dr 


(Al) 


where £. satisfies Jq^^) = 0. Assume that the flux goes to zero at the walls of the 
reactor so that 


(pU v z,t) r (p(r , z, t)J Q (r £.) dr (A2) 

is the finite Hankel transform of the flux. Furthermore, since S6>/ Sr | r _ Q = 0 and 
6(R, z, t) = 0 R (z,t), 


R 

f rp j- (r J 0 (^ r ) dr = aR li e R( z > t ) J i(^i R ) - M) 


To apply a finite Hankel transformation to equation (10), multiply by rJ Q (r and 
integrate from zero to R. Then applying the previous definitions to the result gives 

(— + Vq — \ 0 = Q 0 <p - a 0 + oiR^. 0 R (z, t)J^(^R) 

\St u Sz/ 


(A3) 


The transformed boundary conditions are 


e(%,z,o) = e 0 (£.,z) 
0,t) = 0 


(A4a) 
(A 4b) 
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To solve equation (A3) subject to the boundary conditions (eqs. (A4)), apply a Laplace 
transformation (ref. 22) to equation (A3). Define 


^(%, Z ,P) = f e' pt (p(^,z,t) 


dt 


and 


noting that 


0(%, Z ,P)= f e -pt 0(^,z,t) dt 
<Jf\ 


f e_Pt af dt = z) + p ^ ( %> z,p) 


To apply the transform to equation (A3), multiply by e p ^ and integrate over 1 
zero to infinity; then employing the previous definitions gives 




dz v v o / v 0 v 0 


The solution to this equation is simply 


r 0 


0(^,2, P) = e 


■(P+«|j ) z / v 0 


f Q 0 = 

Jo v 0 


(p+«C)y/v 0 

p)e 1 U dy 


/ Z i _ (p+«?i)y/v 0 f 2 

v 0 Jo 


= (p+a^)y/v 0 

X 0 R (y, p)e dy + 9(^, 0, p) 


a 


(A 5a) 


(A 5b) 


from 


(A6) 


SjRJ^R) 


(A7) 
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Equation (A 4) shows that 0, p) = 0. Equation (A7) may be cast into a slightly different 

form by changing variables twice, first to y/v Q and then to x = (z/Vq) - (y/vg). The 
result is 


Jc\ 


■»z/v n o 

, °_ -(P+aq)x 

0(^,z,P)=Q o / <p(%,z-VoX,p)e dx 


z/v n 2 

-(p+«lf)x 

0 o (4i,z-VoX)e dx 


/ * /v 0 


o^RJ^R) 


/ 


z/ y n 2 

°- -(p+a£. )x 

0 R (z-VgX,p)e 1 dx (A8) 


In order to find 6( r, z,t), ^(| i5 z, p) must be inverted. Now the inverse of the 
Laplace transformation is simply a contour integration over the variable p. Denote this 
integration by the symbol Assuming that the operation may be interchanged 

with the integration over x in equation (A8) it can be seen that 


St 


-lf= 


i ,z,P)] = Q( 


W v ( 


sr 


Z - V ()X, P)e 


■(P+ffq)x 


dx 


Wv„ 


<T 


Ofaz-v&ie 


-(P+a£i)x 


dx 


/ / v 0 

^-1 


^ R (z-V(jX, p)e 


-(P+«^)x 


dx (A9) 


<?(%,z-v qX, P)e 


Now ^ _1 [0(^,z,p)] = e(^,z,t) and ST 1 } 

-a$x L 

xe u(t-x) where u(t-x) is the unit step function defined by 


-(P+a|px 


<p(£it Z-VgXjt-X) 
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u(t-x) = 


fO t < x 

ll t > X 


(A10) 


iL 1 _ -Qf|fx 

? [ 0 o^i» z_v 0 x ^ e e_P J = 0 (A’ Z “ V 

6(t-x) is the Dirac delta function and 


Also note that -S?" 


Qx)e 1 6(t-x)u(t-x) where 


1 0 


<e~ 0 R (z- Vo x,p)e 





t?)x 
=>1 } 


0r(Z- v qX, t- x ) e u(t-x) 


Using the previous relations reveals that the application of the inverse of the Laplace 
transform to equation (A 8) yields 


/ z/v 


°_ -a£x 

z-VqX, t-x)e n u(t-x) dx 


r z/v o 

*/r\ 


+ I V^i’^O 5 ^ 6 

y 0 


a h x 6(t-x)a(t-x) dx 


/ z/v n 2 

0 - 0 !|fx 

0 r (z-VqX, t-x)e u(t-x) dx (All) 


Let t be the time it takes for the fuel to traverse the reactor, that is, r is the transit 

a <* 

time across the reactor. Consider the time interval 0 < t < r a for axial positions in the 
reactor where t < z/vq. Then the unit function u(t-x) is zero for values of x greater 
than t and 0^, z,t) becomes 
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Ill IHIIIIIHI 


f -a^x dx 

e(%,z,t)=Q 0 # <p(%, z-v^t-xje A 


-a I fx 

^f^z-v^e 8(t-x) dx 


/ 

‘'O 

Jc\ 


-a£?x 

0 R (z- Vo x,t-x) e dx 


0 < t < T„ 


t < z/v 


0 


/•to 

Now by definition, / f(x)6(t-x) dx = f(t) provided < t < tg. Thus, the previous 


equation finally becomes 




9 9 

-ap _ -art 

0(l i ,z,t)=Q o 1 Z - V() x, t-x)e dx + 6 0 (4., z-v 0 t)e 


^.RJ^R) / 


0 R (z-VQX,t-x)e dx (A 12) 


Next consider the case when 0 < t < t r and t > z/ V q. In this case the unit function is 
unity over the entire range of the variable x. However, the integration over the delta 
function does not include t in this case so that the integral is zero. Thus 


e(|.,z,t) = Q 


/• 2/v t 

° ./ 

t/n 


_ -«lfx 


/ z/v 0 


-a^x 

0 r (z-VqX, t-x)e dx 


0 < t < r, 
t > z/v f 


a’ 


(A 13) 
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When t > r a , the unit function has value unity over the entire range of x and the integral 
over the delta function never contains t. Therefore ~8 is exactly as in equation (A 13), 
that is, 


r z/v < 

0(^,z,t)=Q o / 


_ -a;£?x 

z-VqX, t-x)e dx 


/ z/v Q 


-a^fx 

e R (z- v QX,t-x) e dx 


t > t (A 14) 


To obtain 0(r, z, t), apply the inversion operation for the finite Hankel transformation, 
that is, 



e(q, z, t) 


J 0<rij> 


( 11 ) 


where the appropriate relation for B( z,t) must be obtained from equations (A12), (A13), 
and (A 14). 
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APPENDIX B 


APPLICATION OF FOURIER-BESSEL EXPANSION TO THE FLUX 


Consider the kinetic equation for neutrons in the form of equation (7), that is, 


^th^ “ ( 6k 0 + + 


From appendix A, it is known that 


-2 


ii/rM 

r 9r \ 9r / 


+ 

9z 2 J 


- y6cp 


0(r, Z,t)= — 6(L, Z,t) 


R 



J tM> 


and 


f z/v ° .«& 

&U i ,z,t)=Q 0 I z-VqX, t-x)e 1 u(t-x)dx 


(Bl) 


( 11 ) 


+ 


/ z/v n 2 

^g^z-v^e 6(t-x)u(t-x) 


dx 


/ z/v 


0 .2 
-«i-x 

0r(z- v qX, t-x)e u(t-x) dx 


Assume that the flux can be expanded in the form <p( r, z,t) = ^ $^(t)Z n (z)(Rj (r) where 

n J 

Z n (z) = sin (e n z), (ft.(r) = J Q (^r), e n = nir/a, £. = za/R, and J Q (^) = 0, Then 
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= jf R r^r^t^r^) dr 

- X /' RrJ 0«( r > r 0«i r ) cir 

m, Z u 

-T ^ *L< t > z m< z >[fi<'';>] 2fi u 

m, l 

- Y E 

4 m 

Using the previous relations in equation (Bl) results in 


*th Z = < 6k 0 + P L> 

n J 


Yj <|, n (t ) Z n (z)(R j (r ) 

n J 


+ p l b " 2 



jz (z) — — 

d*R,(r) 

r J 

d 2 Z (z)l 
+ tfl tr) n > 

r dr 

dr 

-r uij / 

dz 2 , 



-«sfx 

x Z m (z-VQx)e u(t-x) dx + 



_ - a ?ifx 

eu it z-VgX)e 6(t-x)u(t-x)dx 


+ o'^RJjdjR) 



-a^x 

0r(z- v oX, t-x)e u(t-x)dx> 


(B2) 
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Now 




2 2 2 2 2 2 
and d Z n (z)/dz = e*Z n (z). Define B n j = e n + ^ ; then equation (B2) may be written as 


*th X ^i <t)Z n (l)fl J (r) = X [ 6k 0 + P L B ' 2 < b2 - <>] 




n,J 


y 4(t)z n (z) ( R j i 




y ) ^ x (t)Z n (z)(R j (r)(ft i (r)‘jQ 0 

/ J 

n,j 


/ z/v Q 

4 


(t-x) z m (z- v (jX) e 




x u(t-x) dx + 


R 


r^i 

iM 2 i 


z/v. 


-“£i X 

0(^, z-VgX)e 6(t-x)u(t-x) dx 


/ */ v n 
0 

u 


-« A 


0 r (z-VqX, t-x)e 1 u(t-x) dx 


J 


(B3) 


Multiply equation (B3) by Z^ (z) and integrate from zero to a noting that 

£ Z,(z)Z m (z) d z=c m 6; ;m ; 
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l th 2 4|(t)c z (R j (r)= E [5k 0 + P L B- 2 (B 2 - B 2 j)J ^(t)c^(r) 
3 3 


■u 




f z/v 0 

r° 1 4 


n ,J 

m. l 


( t - x )Z m ( Z -v 0 x) 


-a£x 2 

x e u(t-x) dx dz + 


r2 ^iM 2 4 


j z/v ( 


0^, z_v 0 x ) e 


•a^fx 


x 6(t-x)u(t-x) dx + 


/ z / v n 2 

-a^fx 

SpCz-Vf^x, t-x)e 1 u(t-x) dx 


Now if the z and x integrations are interchanged, the result is 

rSL /-z/v n 1 /*a r z 

/ dz / u f(z)g(z-v n x) dx = — / dz / f(z)g(z-s) ds 
* / 0 * / 0 ^ v Q * / 0 y 0 


(B4) 


= — f ds /* f(z)g(z-s) dz 
V n * / 0 “o 


r a/v n r 2i 

/ 0 dx f i{z)%{z-x tf.) dz 

•'O ‘'vqX u 


= f a dx f f(z)g(z-v n x) dz 

* / 0 *v qX ^ 


where f and g are arbitrary functions and s = VqX. Define 


.a 


I W v 0 x > =-r L, V z > Z l( z ) Z m (s -Y) dz 

c z v tr 


(B5) 
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Then equation (B4) becomes 


*th £ij(t)«yr)= E[6k 0+ P L B- 2 (B 2 - B 2 .)J (t)fftj(r) 
3 3 


- r y 4 (t) J dx(R.(r)(R i (r)<Q 0 -fj n (t 

n, j 0 

m, i 


-a%fx 

t_x)K nZm^ v O x ) e “(*-*) 


c l R " 


f - 

/ ^O^i’ z_v o x ) e 

m 2 4 


Z i (z)Z n (z)6(t-x)u(t-x) dz 


2 ^i f 

c i w M 1^ 

v Cr 


«lfx 


^R^ Z ~ V 0 X ’ t-x)Z^ (z)Z n (z)e 1 u(t-x) dz 


(B6) 


Next, multiply equation (B-6) by r& k (r) and integrate from zero to R noting that 


/ R 2 

r W w o^ r > dr =y[ J i("i>] \ 


The result is 
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Y fiM 2 = h + p l b ' 2 < b2 - >&>] 


- y 


E $^(t) / dx / r(fyr)<ft (r)(R.(r) 

J 0 Jo ’ 

n, j 

m, i 


M Q 0*L< t - x > I W' , 0 !t,e aSiXu(t - x) 


c z R 


£ 


L_ [ 

l^ij] 2 gX 


Z z (z) z n (z)0g(4 i , z - V o x) e 6(t-x)u(t-x) dz 


/ a 

v o x 


0 R (z-VgX,t-x)Z z (z)Z n (z)e u(t-x) dz 


(B7) 


Defining 


M ijk- 2 


2 />R 

Jo 


Ey^J] 


gives the desired equation for the coefficients 


r (Rjj (r )(R.j (r )(R i (r ) dr 
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L <m { 

n,j 




& * «<*-*) 


V_ 


c z R 


w 4 


_ -a£ 2 x 

Z^(z)Z n (z)0o(^, z-VQx) e 5(t-x)u(t-x) dz 




2 <^i f 

^ 4 


-a^x 


0 R (z- v o x ,t- x )Z z ( z ) z n (z)e 1 u(t-x) dz> (B8) 




APPENDIX C 


VALIDITY OF ONE-TERM APPROXIMATION 

To investigate the one -term approximation, the steady -state flux amplitude and the 
corresponding temperature rise at r = 0, z = a (exit temperature) will be calculated. IE 
the thermal conductivity is not too large, the axial and radial fluxes can be considered as 
independent quantities. Two cases will be investigated. First, the thermal conductivity 
will be set equal to zero and the exact solution to the flux and temperature equations will 
be obtained. The conditions under which sin e^z is an appropriate description of the axial 
flux will be found from this case. Second, assuming a single axial term is appropriate, 
a two -term radial approximation will be used in the coupled equations and the conditions 
under which the two-term flux and exit temperature agree with the one-term solutions 
will be investigated. 

The basic equations are the following steady-state flux and temperature equations: 



1 

r 



de 

3z 



1 

r 



(Cl) 


(C2) 


Case 1: ot = 0 


In this case, it is assumed that the radial temperature follows the flux to obtain 
0(r, z) = (fl^(r)u(z) and cp(r, z) = <%( r )<P s ( z ) where (ftj^r) = JQd-^r). Noting that 

B 2 = e 2 + | 2 and that 


I A 

r dr 



= - 


means equations (Cl) and (C2) may be written as 



+ (5k Q - = 0 
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v 0 du 
Q 0 dz 


tfl 1 = ^ g (R 1 


Multiplying these equations by r(R 1 (r) and integrating from zero to R yield (Note that 

-L 

<r 2 / 2 )[j 1 (, 1 )] 2 = r n r(R^(r)(Rj(r) dr and = J rdl^(r)(Rj(r)(Jlj(r) dr 

x 2/R 2 jj 1 (z' 1 )j 2 . ) 


-Pt /d 2 ^ s 2 . 

— ( ;~ + e l (p s] + 6k 0^s ‘ Y V( Ps M 1U ~ 0 

B Z \dz 2 


(C3) 


and 


Q 0 dz 


(C4) 


Using equation (C4) in equation (C3) gives 


d yA t A-!>Li = o 

dzl dz 2 2 


(C5) 


2 2 


where \i s P L /B , y ^ = yM.^, and Pj = 6k Q + h e^. Integrating equation (C5) gives 


.2 d 2 u 2 r l 2 , . . 

h — — + p?u — - u = f = constant 

dz 2 2 

p p 

Let y = du/dz and y dy/du = d ii/dz ; then the previous equation becomes 


l 2 dy 2 r l 2 , 

h y ^ + p^ u = f, 

du 1 2 1 


Integrating this yields 


2 2 P 2 u 2 y^ 3 
h 2 L_ = - J_ + _L_ + f u + f 

ri o ‘ 
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where fg is a constant. Using the definition of y results in 



2 2 c 

V /1 U 


6 


+ f l u +f 2 


(C6) 


Now the axial temperature rise is a minimum and is equal to zero at the inlet and is a 
maximum at the outlet, so that du/dz = 0 when z = 0, a and u = 0 when z = 0. Thus, 
fg = 0 and 


f 


1 " 


p?u 
^1 ex 

2 



6 


where vi = u when z = a. Thus, equation (C6) finally becomes 

CA 



p 

If Tj/T^ ex = sin Xp the previous equation can be cast into the form 


-sin’ yTj/T^ 

T -o / 3 j 

I *1 

V 3 - T lex J 

^1 - K 2 sin 2 x 1 

4 

1 

where K 2 = T Jex /(3 - T lex ). When r . 

r l = T lex’ * = * 
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Equation (C7) contains a complete elliptic integral of the second kind; an approximate 
value of the integral can be obtained if K 2 « 1. Then 


^1 - K 2 sin 2 


- 1/2 


=* 1 + 


v 2 .2 v 

K sin x 1 


and 


a = 7 r 


V 3 - T lex 



(C8) 


o 

Equation (C8) must be solved for T lgx when K is small; the solution is 

T lex ==; 26kQ/h 2 e^ where it has also been assumed that 6k Q « h 2 e^. Using the defin- 
ition of T^ ex results in the steady-state exit temperature 


u 


ex 


28k 0 

>1 


(C9) 


2 

Note that two things have been assumed: first, that K « 1 and, second, that 
6k q « h 2 e 2 . These assumptions are quite good' in most cases of practical interest. The 
first assumption holds provided that u « 1 which is obviously true if ok Q « 1. 
This is essentially the same as the second assumption. 

Now the one-term steady-state equation for a = 0 and t > r a is simply (see 
eq. (17)) 


6k, 


*0 = 


K i ^ aK m( v o x ) dx 


(CIO) 
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where K^=yQ 0 M^^^ and J' a Kill^o’O dx = r^/n. The steady- state temperature is 
given in appendix D as 


<V r , z) = - 


Wo' 1 *!) L2 


k ) 2 * <y 


a 


a sin e-jZ — — l cos e^z - e 


■ a h z/v o\ 


which, when evaluated at r = 0, z = a and using equation (CIO) for 4 ?q yields 

26k, 


0 n (O,a) = u ex = 


"0 


(D7) 


(Cll) 


y 1 


The equality of equations (C9) and (Cll) indicates that a one-term axial approximation is 

2 2 

adequate to describe the system provided that 6k Q « P h q/B ■ 

Case 2: a + 0 


The first two radial term equations for the steady-state flux amplitude are 

r g ™ 

0 = $}6k 0 - Q 0 y ($}) M 111 g} 11 + $! M 121 G 111$1 + $ 1 $ 1 M 211 G 111 + $ 1 M 221 G 111 $ 1 

(C 12) 


0 = 


( b2 - B L) + 6k ol - %yH M 112 G hA + *1*? M 122 G 111 

+ ^i^i M 212 G 111 + ^l M 222 G lll^l) ( C13 ^ 


where 




rtfl^rjtft^rj^r) dr 


and 
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g: 


nZm 


r 

J(\ 


K 


nZm 


(VqX)! 


~ a & 


dx 


Let 


F 12=-^( b2 - B 12) + 51I 0 


multiply equations (C12) and (C13) by 4>j/ and define /3 = Then 

6k 


,2/U 


yQ 


- = 4 i M m G m + * 

0 


1 3M 112( G 111 + G lll) + * 1 ? 


2m 112 G U1 


f> — = *1 M U2 G 111 + 4 


Q 


0 


l^mKn ♦ G iii) + *}(> 


2-iwr p 2 

m 222 u 111 


(C 14a) 


(C14b) 


Note that, when 0 = 0, equation (C14a) yields the one-term approximation for that is, 

i 8k « 




1, 1-term 


(C15) 


^0 M 111 G 111 

Solving equations (C14) for 4>J and equating the results yield the following equation: 

^ 3 ( M 122 G 111 F 1 2 ) + ^mHll + G Ul) F 12 ‘ M 222 G m 6k Q] 


^INhx g x 


111 F 12 " M 122 


v-1 


( G lll + G lll) 6k J - M 112 G lll 6k 0 = 0 < C16 > 


where the common factor (tQq) has been multiplied out. This cubic equation must 
now be solved for j3. It can be shown by a straightforward calculation that, if 

i 2 2 P i 2 

5k 0 « |P L (q- q)/B‘ *\ and « 1, then there is only one real root for 0 which 

is 


0 = 0. 0005 1 + 2 * 
3 


(C17) 
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where = r a a (^\ - ^)/ ^ « 1. Using equation (C17) in equation (C14a) gives 


1, 2 -term 


yQ^lllCin 1 + M 112 L | G lll 


M m V g} u 


where terms in have been neglected. Now 




and M 112 /M 1;11 = 0. 185 (l/^) where = J 1 (^ 1 )/J 1 (t' 2 ). Thus 


*1. 2-term = 0 - 9998 


9 -_ = 0.9998 $} 1 _ term 

^ Q 0 M 111 G 111 


2 1 

and ^ 0. 00049 <f>^ i_t erm - Thus, the one-term approximation is accurate to within 
0. 05 percent for <$. ’ 
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APPENDIX D 


INITIAL AND STEADY -STATE FLUX AND TEMPERATURE DISTRIBUTIONS 

The initial flux may be obtained from equation (7) with dcp/ 3t = 0 and 6k = 0. Then 


I±(r!f?i) + if® + B*,. 0 =0 


r 3r \ 3r 


(Dl) 


3z 


n 

where B^ = - S a )/D is the physical buckling. The solution to the equation is 


<P 0 = # 0 sin ( e i z ) J o^i r ) 


P 2) 


which can be demonstrated by using equation (D2) in equation pi). Here, <3 ?q is an 
amplitude which depends on the particular initial conditions in the reactor. Furthermore, 
B 2 = B 2 = e 2 + | 2 for this solution to hold. 

The initial temperature distribution can be obtained from equation (10) with 30 / 3t = 0, 
that is, 


!V% 0+ jl 

3z Vq v 0 r 9r \ 9r / 


P3) 


The boundary conditions are 0Q(r, 0) = 0, 0q(R, z) = 0, and dO^/dr | r _Q = 0. Then applying 
a finite Hankel transform (see appendix A) to equation (D3) results in 


<u 0 _ Q 0 4. 0 R 2 ]^) 


dz 


2 

6., sin e.z i 0 Q 

v 0 


P4) 


where 0q(^, z) is the transformed temperature and equation (D2) is used for the flux. 
Since 0 q(£., 0) = 0 from the first boundary condition, the solution of equation (D4) is 

%*o r2 |j 
v o 



0 o u v z ) = e 


■a^z/v 


0 


/ 
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Performing the integration in equation (D5) gives 






\a(f) + <rA a ) 2 


a if sin £ jz - — (cos tjZ - e '“ 4 i z/v o\ 




Using the inversion formula finally yields 


R 2 L-J [^(-i)]' 


(D6) 


where 


0 o (r, z) = 


QqVo^i) 

H f + WrJ 


at* sin €l z - — (cos €l z - 


(D7) 
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APPENDIX E 


EVALUATION OF STABILITY CRITERION 


The following integral will be investigated: 


I = — /*°° S(x) sin fix dx 

7T *A) 


1 -a£?x 

where S(x) = G(x)/G(0) and G(x) = y M ni(P c p)" Kj^ 1 (vQx)e u(T a -x). Integrating by 
parts and using the fact that G(<») = 0 results in 


2fi /•» 


ttG(O) */0 


f G(x) 

Jo 


cos fix dx 


Using the definition of G(x) gives 


r 


-a^x 

K m ( v oX)e cos fix dx 


Now when (xHt « 1, the exponential term adds very little to the integral. For con- 

J. cl 

venience, therefore, the exponential is excluded and 


I « — f a K 
4 y 0 


IH^qX) 


cos fix dx 


K , . . (v„x) = i (l + - cos — + I cos — 

"V 3 T a 3 T a, 


so that equation (E4) becomes 


1=1 


sin fir„ 


'if . nfef - n 2 


I 



It follows immediately that I >: 0 if S2 < Sir / t . For values of 52> Sir / t , I oscillates 

a d. 

with sin fiT . Equation (E 5) specifically precludes O = n / r or S2 = 2ir/r ; when 

a a a 

either of these values are used in equation (E4), the value of I is greater than zero. 
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